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ABSTRACT 

Due to the chaotic nature of the Solar System, the question of its dynamic long-term stability 
can only be answered in a statistical sense, for instance, based on numerical ensemble integrations 
of nearby orbits. Destabilization of the inner planets, including catastrophic encounters and/or 
collisions involving the Earth, has been suggested to be initiated through a large increase in 
Mercury’s eccentricity (ex), with an estimated probability of ^1%. However, it has recently been 
shown that the statistics of numerical Solar System integrations are sensitive to the accuracy and 
type of numerical algorithm. Here I report results from computationally demanding ensemble 
integrations {N = 1,600 with slightly different initial conditions) at unprecedented accuracy 
based on the full equations of motion of the eight planets and Pluto over 5 Gyr, including 
contributions from general relativity. The standard symplectic algorithm used for long-term 
integrations produced spurious results for highly eccentric orbits and during close encounters, 
which were hence integrated with a suitable Bulirsch-Stoer algorithm, specifically designed for 
these situations. The present study yields odds for a large increase in Mercury’s eccentricity 
that are less than previous estimates. Strikingly, in two solutions Mercury continued on highly 
eccentric orbits (after reaching ex values >0.93) for 80-100 Myr before colliding with Venus 
or the Sun. Most importantly, none of the 1,600 solutions led to a close encounter involving 
the Earth or a destabilization of Earth’s orbit in the future. 1 conclude that Earth’s orbit is 
dynamically highly stable for billions of years in the future, despite the chaotic behavior of the 
Solar System. 

Subject headings: celestial mechanics — methods: numerical — methods: statistical — planets and 
satellites: dynamical evolution and stability 


1. Introduction 

One of the oldest and still active research areas 
in celestial mechanics is the question of the long¬ 
term dynamic stability of the Solar System. After 
centuries of analytical work led by Newton, La¬ 
grange, Laplace, Poincare, Kolmogorov, Arnold, 
Moser etc. (Laskar 2013), the field has recently 
experienced a renaissance due to advances in nu¬ 
merical algorithms and computer speed. Only 


since the 1990s have researchers been able to in¬ 
tegrate the full equations of motion of the Solar 
System over time scales approaching its lifetime 
(± ^5 Gyr) (Wisdom & Holman 1991; Quinn et al. 
1991; Sussman & Wisdom 1992; Saha & Tremaine 
1992; Murray & Holman 1999; Ito & Tanikawa 
2002; Varadi et al. 2003; Batygin & Laughlin 2008; 
Laskar & Gastineau 2009; Zeebe 2015), and ex¬ 
ceeding Earth’s future habitability of perhaps an- 
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other 1-3 Gyr (Schroder & Smith 2008; Rushby 
et al. 2013). Looking ahead, detailed exoplanet 
observations (e.g. Oppenheimer et al. 2013) will 
enable similar integrations of planetary systems 
beyond our own solar system. 

Importantly, in addition to CPU speed, a sta¬ 
tistical approach is necessary due to the chaotic 
behavior of the Solar System, i.e. the sensitivity of 
orbital solutions to initial conditions (Laskar 1989; 
Sussman & Wisdom 1992; Murray & Holman 
1999; Richter 2001; Varadi et al. 2003; Batygin & 
Laughlin 2008; Laskar & Gastineau 2009; Zeebe 
2015). To obtain adequate statistics, ensemble 
integrations are required (Laskar & Gastineau 
2009; Zeebe 2015), that is, simultaneous inte¬ 
gration of a large number of nearby orbits, e.g. 
utilizing massive parallel computing. Chaos in 
the Solar System means small differences in ini¬ 
tial conditions grow exponentially, with a time 
constant (Lyapunov time, e.g. Morbidelli (2002)) 
for the inner planets of ^^3-5 Myr estimated nu¬ 
merically (Laskar 1989; Varadi et al. 2003; Baty¬ 
gin & Laughlin 2008; Zeebe 2015) and ~1.4 Myr 
analytically (Batygin et al. 2015). For exam¬ 
ple, a difference in initial position of 1 cm grows 
to -1 AU (= 1.496x10“ m) after 90-150 Myr, 
which makes it fundamentally impossible to pre¬ 
dict the evolution of planetary orbits accurately 
beyond '^lOO Myr (Laskar 1989). Thus, rather 
than searching for a single deterministic solution 
(conclusively describing the Solar System’s future 
evolution, see Laplace’s demon, Laplace (1951)), 
the stability question must be answered in prob¬ 
abilistic terms, e.g. by studying the behavior of a 
large number of physically possible solutions. 

Until present, only a single study is available 
that integrated a large number of Solar System 
solutions based on the full equations of motion 
and including contributions from general relativ¬ 
ity (Laskar & Gastineau 2009). Using a sym- 
plectic integrator throughout, Laskar & Gastineau 
(2009) integrated 2,501 orbits over 5 Gyr and 
found that Mercury’s orbit achieved large eccen¬ 
tricities (>0.9) in about 1% of the solutions. Fur¬ 
thermore, they suggested the possibility of a col¬ 
lision between the Earth and Venus. Given that 
only one study of this kind exists to date, sev¬ 
eral questions remain. For instance, are the nu¬ 
merical results and statistics obtained sensitive 
to the numerical algorithm used in the integra¬ 


tions? Furthermore, is the symplectic algorithm 
used throughout a reliable integrator for highly 
eccentric orbits and during close encounters? If 
not, what are the consequences for Solar System 
stability? In particular, considering long-term sta¬ 
bility and planetary habitability, what are the con¬ 
sequences for the evolution of Earth’s future orbit 
over billions of years? Addressing these questions 
is the focus of the present study. 

2. Methods 

I have integrated 1,600 solutions with slightly 
different initial conditions for Mercury’s posi¬ 
tion based on the full equations of motion of 
the eight planets and Pluto over 5 Gyr into the 
future using the numerical integrator packages 
HNBody (Rauch & Hamilton 2002) and mercury6 
(Chambers 1999). Relativistic corrections (Ein¬ 
stein 1916) are critical (Laskar & Gastineau 2009; 
Zeebe 2015) and were available in HNBody but not 
in mercuryG. Post-Newtonian corrections (Soffel 
1989; Poisson & Will 2014) were therefore imple¬ 
mented before using mercuryS (see Appendix). 
Thus, all simulations presented here include con¬ 
tributions from general relativity (GR). To allow 
comparison with previous studies, higher-order ef¬ 
fects such as asteroids (Ito & Tanikawa 2002; Baty¬ 
gin & Laughlin 2008), perturbations from passing 
stars, and solar mass loss (Ito & Tanikawa 2002; 
Batygin & Laughlin 2008; Laskar & Gastineau 
2009) were also not included here. 

The initial 5-Gyr integrations of the 1,600- 
member ensemble were performed with HNBody 
(Rauch & Hamilton 2002) using a 4th-order sym¬ 
plectic integrator plus corrector with constant 4- 
day timestep (At, Table I). The ensemble com¬ 
putations required '^6 weeks uninterrupted wall- 
clock time on a Cray CS300 (~1.7 million core 
hours total), plus up to ^4 months for individ¬ 
ual runs at reduced timestep (see below). For the 
symplectic integrations, Jacobi coordinates (Wis¬ 
dom & Holman 1991) were employed rather than 
heliocentric coordinates, as the latter may under¬ 
estimate the odds for destabilization of Mercury’s 
orbit (Zeebe 2015). 

All simulations started from the same set of 
initial conditions (Table 2), except Mercury’s ini¬ 
tial radial distance was offset by 1.75 mm between 
every two adjacent orbits. The largest overall off- 
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Table 1: Summary of the 5-Gyr simulations. 


bm “ 

Algorithm'’ 

At 

(days) 

N 

^Collisions'’ 
M-V MS 

< 1.00 

4th-sympl. 

4 

1600 


> 0.55 

4th-sympl. 

1 

28 


> 0.70'^ 

4th-sympl. 

1/4 

10 


> O.SO'’ 

BS 

adaptive 

10 

7 3 


°‘€.M = Mercury’s eccentricity. 

^4th-sympl. = 4-th order symplectic (HNBody). BS = 
Bulirsch-Stoer (mercury6). 

'’’At-V; Mercury-Venus, M-S\ Mercury-Sun. 
"^Switched to 1/4 day-symplectic in case of large 
\AE/E\ variations at 1 day. 

'’Switched to BS in case of ajvi shifts and/or large 
\AE/E\ variations at 1/4 day-symplectic. 

set was 1,599 x 1.75 mm ~ 2.80 m, well within 
the uncertainty of our current knowledge of the 
Solar System. Initial conditions for all bodies 
in the 5-Gyr runs (before offsetting Mercury) 
were generated from DE431 (naif . jpl .nasa. 
gov/pub/naif/generic_kernels/spk/plELnets) 
at JD2451544.5 (01 Jan 2000) using the SPIGE 
toolkit for Matlab (naif.jpl.nasa.gov/naif/ 
toolkit.html) (Table 2). The small offsets in 
Mercury’s initial position randomized the initial 
conditions and led to complete divergence of tra¬ 
jectories after ~100 Myr (see Eig. 4). 

In case Mercury’s eccentricity increased above 
certain threshold values during the simulation (Ta¬ 
ble 1), At was reduced but held constant after that 
until the next threshold was crossed, if applicable. 
For example, a 4-fold reduction in stepsize, twice 
throughout the simulation (i.e. 4 —>• 1 —>• 1/4 day) 
typically kept the maximum relative error in en¬ 
ergy \AE/E\ = \{E{t) — Eo)/Eo\ and angular mo¬ 
mentum (|AL/L|) of the symplectic integrator be¬ 
low 10“^° and 10“^^ (Fig. 1). These steps in 
At corresponded to cm thresholds of about 0.55 
and 0.70 (Table 1). All results of the symplectic 
HNBody integrations were inspected and (if appli¬ 
cable) restarted manually with smaller At at the 
appropriate integration time using saved coordi¬ 
nates from the run with larger At (automatic step- 
size reduction was not possible because HNBody’s 
source code is not available). 

Importantly, high cm solutions associated with 


qm shifts and/or large \AE/E\ variations (usu¬ 
ally during close encounters) were integrated using 
the Bulirsch-Stoer (BS) algorithm with adaptive 
stepsize control of mercuryS, including GR contri¬ 
butions (see Appendix) and specifically designed 
for these tasks (Chambers 1999) (the symplectic 
algorithm produced spurious results in these sit¬ 
uations, see Figs. 2, 6). Specifically, the sym¬ 
plectic HNBody integrations (1/4 day) were in¬ 
spected for aM shifts and/or large \AE/E\ vari¬ 
ations and restarted with mercuryS’s BS algo¬ 
rithm at the appropriate integration time using 
saved coordinates from the symplectic HNBody run. 
The only comparable study to date (Laskar & 
Gastineau 2009) integrated 2,501 orbits using a 
symplectic algorithm and an initial At ~ 9 days. 
While Laskar & Gastineau (2009) also reduced the 
timestep depending on cm (though with gener¬ 
ally larger maximum errors, see below), high cm 
solutions and close encounters were integrated us¬ 
ing the symplectic algorithm throughout (Laskar 
& Gastineau 2009). However, for highly eccentric 
orbits the sympletic method can become unsta¬ 
ble and may introduce artificial chaos, unless At 
is small enough to always resolve periapse (Rauch 
& Holman 1999). In the Solar System, At must 
hence resolve Mercury’s periapse with the highest 
perihelion velocity (vp) among the planets (and 
increasing with e, as Wp oc (1 -I- e)/(l — e)). 

3. Results 

The vast majority of solutions obtained here 
showed a stable evolution of the Solar System 
and moderate cm values (eA<<0.6 in 99.3% of all 
runs). In 10 out of 1,600 solutions (^0.6%), Mer¬ 
cury’s eccentricity increased beyond 0.8 (see be¬ 
low). Importantly, the maximum relative error in 
energy \AE/E\ and angular momentum (|AL/L|) 
of the symplectic integrator was typically held be¬ 
low 10“^° and 10“^^ even at ~ 0.8 by e.g. a 
4-fold reduction in stepsize twice throughout the 
simulation (Fig. 1). Gritically, once cm increased 
beyond ^0.8 (resulting in close encounters with 
Venus and shifts in Mercury’s semi-major axis, 
awi), the Bulirsch-Stoer integrator of mercuryG 
was employed. For example, in run ^0299 (R0299 
for short) close encounters between Mercury and 
Venus occurred at t cs 4.3042 Gyr at which time 
they approach their mutual Hill radius (Chambers 
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HNBody, 4th-order sympi, Run #0649 



Fig. 1.— Example of a symplectic multi-billion year 
integration of the Solar System with HNBody (run 
#0649). Results shown are from symplectic integra¬ 
tor only (4-th order with corrector and Jacobi coordi¬ 
nates), not from Bulirsch-Stoer integrator (see text). 

(a) Relative energy error, |AJ?/J?| = \{E{t) — Eo)/Eo\. 
After the two reductions in integrator timestep (ti, 
i = 1,2), solid and dashed lines show \AE/E\ rel¬ 
ative to E{ti) and relative to absolute Eo{to = 0), 
respectively, (b) Relative angular momentum error, 
|AL/L| = \{L{t) — Lo)/Lo\. (c) Change in Mercury’s 
and Earth’s semi-major axes, Aa/a — {a{t) — ao)/ao. 

(d) Mercury’s and Earth’s eccentricity {cm, eg), (e) 
Inclination. Note that \AE/E\ < even at 

Cm — 0.8. When oscillations in \AE/E\ and/or shifts 
in aM occurred (cm > 0.8), integrations were con¬ 
tinued with mercury6’s Bulirsch-Stoer algorithm (see 
Figs. 2, 5, 6). 

et al. 1996); 

0.0053 AU, (1) 

2 \ 3mi / 


fas (Myr) 


0 0.5 1 1.5 2 2.5 3 3.5 



4.304 4.3045 4.305 4.3055 4.306 4.3065 4.307 
Time (Gyr) 


Fig. 2.— Symplectic and Bulirsch-Stoer (BS) inte¬ 
gration of run #0299 at high cm- The BS integra¬ 
tion (mercury6, green lines, approximate relative x 
and V error per step = 10“^®) was initiated at tes = 0 
using the coordinates of the symplectic run (HNB = 
HNBody, purple lines, 4-th order, timestep = 1/4 day), 
(a) Relative energy error, \AE/E\. (b) and (c) Mer¬ 
cury’s semi-major axis, aj^ (note different y-scales). 
(d) Mercury’s eccentricity, cm- Note spurious drop in 
aM at tss — 340 kyr in the symplectic integration 
(arrow), which ultimately causes rapid destabilization 
of Mercury’s orbit (see text). In contrast, aM is es¬ 
sentially stable in the BS integration, which properly 
resolves Mercury’s periapse and close encounters as a 
result of adaptive step size control. 

where a’s and m’s are the planetary semi-major 
axes and masses; mi is the solar mass. 

More precisely, close encounters with dmin = 
0.0080 AU ensued a.t Ibs — 340 kyr, where Ibs 
is measured from the start of the BS integration 
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Fig. 3.— Illustration of trajectory separation as Mer¬ 
cury’s orbit approaches instability, (a) cm from two 
symplectic integrations with HNBody of Run 0043 of 
one fiducial (solid blue line, Nr. (1)) and one shadow 
orbit (dashed green line, Nr. (2)) offset by 1.75 mm 
in Mercury’s i-coordinate at t' = 0. (b) Mercury’s 
semi-major axis, (c) Difference in cm (Acai) between 
(1) and (2). At t' ~ 2.2 Myr, close encounters between 
Mercury and Venus ensue (see jumps in um and Aca^ , 
arrows). 

and dmin is the minimum distance between the 
two bodies (Fig. 2). This point in time corre¬ 
sponds to a significant drop of aM to ~0.3869 AU 
in the symplectic integration, while aM in the BS 
integration remains stable. Note that the corre¬ 
sponding symplectic \AE/E\ < 10“® at 340 kyr 
may be misinterpreted to reflect accurate orbit in¬ 
tegration (in fact, if the symplectic At was re¬ 
duced right afterwards, spurious behavior would 
likely go unnoticed). However, this is clearly not 
the case as shown by comparison with BS, which, 
due to adaptive step size control, properly resolves 
Mercury’s periapse and close encounters. Hence 
the relative energy error (say < 10Laskar & 
Gastineau (2009)) is not a sufficient criterion to 
ensure accurate steps in symplectic integrations 
with highly eccentric orbits and close encounters. 


The drop in aM at Ibs — 340 kyr in the sym¬ 
plectic R0299 integration preconditions the sys¬ 
tem for further eM increase, subsequently causing 
a large aM rise and rapid destabilization of Mer¬ 
cury’s orbit (Fig. 2). Specifically, the aM drop is 
followed by an immediate increase in the average 
Cm in the symplectic integration relative to BS, 
which, in turn degrades Mercury’s perihelion reso¬ 
lution (at constant symplectic At) and leads to os¬ 
cillations in \AE/E\ and to further, larger swings 
in aM- In contrast, while in the BS integration 
aM also oscillates, it remains overall fairly stable 
over several million years. Eventually, Mercury 
collides with Venus in BS-R0299 at tss — 15 Myr 
(Fig. 5). Thus, while ultimately both the sym¬ 
plectic and BS integration of R0299 spell disas¬ 
ter for Mercury’s orbit, it involves vastly different 
trajectories and time scales. Notably, Mercury’s 
extended lifetime in the BS integration relative to 
the symplectic algorithm is a common feature (see 
Figs. 5, 6). One reason for this is of course a too 
large (and constant) timestep in the current sym¬ 
plectic integrations. The crux, however, is that 
even if the symplectic timestep was reduced be¬ 
fore |AE/E| becomes too large, spurious results 
could easily be overlooked. 

Note that the spurious aM shift in the sym¬ 
plectic integration of R0299 occurs after only a 
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Fig. 4.— Mercury’s maximum eccentricity per 1 Myr 
bin from the 1,600 symplectic 5-Gyr integrations. Re¬ 
sults shown are from HNBody’s symplectic integrator 
only (Rauch & Hamilton 2002). In ten solutions, cm 
crossed 0.8, after which the integration was continued 
with mercury6’s Bulirsch-Stoer algorithm (Chambers 
1999) (Fig. 5). Note that in R0840 (arrow), cm < 
0.8, no shifts in um occurred, and max\AE/E\ ~ 
8x10"“. 

few 100 kyr during close encounters {tss = 0 
means same coordinates in symplectic and BS in¬ 
tegration) and is therefore not due to intrinsic 
chaos, whose characteristic time scale causes tra¬ 
jectory separation over millions of years. The 
million-year separation time still holds even as 
Mercury’s orbit approaches the regime of insta¬ 
bility in the solutions studied here {cm > 0.8). 
This behavior can be illustrated by following two 
nearby trajectories initiated at high cm (Fig. 3). 
One fiducial and one shadow orbit of R0043 off¬ 
set by 1.75 mm in Mercury’s a;-coordinate were 
integrated using HNBody’s symplectic algorithm. 
Polynomial growth in Acm (difference in cm be¬ 
tween the two orbits) governs Ae^vi’s increase over 
the hrst ~2.2 Myr (cf. Zeebe 2015). However, at 
t' ~ 2.2 Myr, close encounters between Mercury 
and Venus ensue (as in R0299), causing jumps in 
gm and Acm (arrows, Fig. 3). Additional shadow 
runs initiated at t' = —26 Myr and —11 Myr (not 
shown) exhibit no close encounters during run-up 
to t' = 2.2 Myr and indicate exponential trajec¬ 
tory divergence due to chaos after > 8 Myr. Thus, 
the abrupt divergence of trajectories at 2.2 Myr is 
due to close encounters, not intrinsic chaos. 


3.1. High-e^ solutions 

In ten solutions, ex increased beyond 0.8 of 
which seven (three) resulted in collisions Mercury- 
Venus (Mercury-Sun) (Table 1, Fig. 4). Thus, the 
odds for a large increase in Mercury’s eccentricity 
found here (0.6%) are less than previous estimates 
of ~1% (Laskar & Gastineau 2009). Strictly, the 
latter odds were actually 20/2501 = 0.8%. Us¬ 
ing different statistical methods (e.g. Agresti & 
Coull 1998), the 95% confidence interval for the 
two results (0.6% and 0.8%) may be estimated as 
±0.42% (V = 1,600) and ±0.36% [N = 2,501). 
Thus at the 95% conhdence level, one can con¬ 
clude that the two results are not statistically dif¬ 
ferent from one another, which could be remedied 
by performing ensemble runs with even larger N. 
It is important to realize, however, that in or¬ 
der to reduce the 95% confidence interval for the 
same distributions to, say ±0.1% would require 
N > 20,000. Note that the 1,600 solutions used 
initial conditions that differed only by 1.75 mm 
in Mercury’s initial radial distance between every 
two adjacent orbits. Yet the small offsets in Mer¬ 
cury’s initial position randomized the initial condi¬ 
tions and led to complete divergence of trajectories 
after ^100 Myr (Fig. 4). 
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Fig. 5.— 120 Myr of Bulirsch-Stoer integration of ten solutions, initiated when oscillations in |Ai?/i5| and/or shifts 
in aM occurred in the symplectic integrations {cm ^ 0.8, see text), (a) Relative error in energy and (b) angular 
momentum, (c) Mercury’s semi-major axis and (d) eccentricity. Note solutions R0043 and R1231. 
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Fig. 6.— Comparison of symplectic vs. BS integration of all ten high-eAr solutions during the initial phase (cf. 
Fig. 2). In each panel, the top two graphs show log \AE/E\ (left axes), the bottom two graphs show um (right axes), 
(a) Column of five runs with ’initial’ symplectic \/S.E/E\ < 10~® (see text). Note symplectic \l\E/E\ oscillations 
< 10~® (arrows, corresponding to high cm periods) after which symplectic and BS trajectories diverge (indicated by 
different subsequent aM evolutions), (b) Five runs with ’initial’ symplectic \/S.E/E\ > 10~® (see text). 
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The results of the symplectic algorithm were 
used only until shifts in gm and/or large \AE/E\ 
variations occurred in the 1/4 day-symplectic in¬ 
tegration {cm > 0.8, Fig. 4). Subsequently, the in¬ 
tegration was continued with mercuryS’s Bulirsch- 
Stoer (BS) algorithm (Fig. 5) (all integrations 
included GR contributions, see Appendix). As 
is well known, long-term conservation of energy 
and angular momentum in BS integrations (here: 
\AE/E\ < 2x10-9, |AL/L| < 5x10-1°) is gener¬ 
ally worse than in symplectic integrations. How¬ 
ever, as mentioned above, for highly eccentric or¬ 
bits the sympletic method can become unstable 
and may introduce artificial chaos (Rauch & Hol¬ 
man 1999). Thus, in the present case (critical time 
interval < 100 Myr), it was imperative to inte¬ 
grate high Cm solutions and close encounters with 
mercuryS’s BS algorithm (adaptive stepsize con¬ 
trol), specifically designed for these tasks (Cham¬ 
bers 1999). The symplectic algorithm produced 
spurious results in these situations (Figs. 2, 6). 

All ten solutions with ex > 0.8 integrated 
with Bulirsch-Stoer ultimately resulted in colli¬ 
sions involving Mercury; Mercury-Venus: R0299, 
R0397, R0649, R1231, R1320, R1461, and R1594; 
Mercury-Sun: R0043, R0728, and R1322. Strik¬ 
ingly, in two solutions (R0043 and R1231), Mer¬ 
cury continued on highly eccentric orbits (after 
reaching ex values >0.93) for 80-100 Myr before 
colliding with Venus or the Sun (Fig. 5). This sug¬ 
gests the potential existence of non-catastrophic 
trajectories, despite Mercury achieving large ec¬ 
centricities (is there a chance of recovery, evad¬ 
ing detrimental consequences for the Solar Sys¬ 
tem?). In contrast, once ex crossed 0.90-0.95 in 
the symplectic integrations {At = 1/4 day), large 
shifts in ax and/or energy occurred rapidly, of¬ 
ten leading to complete destabilization of Mer¬ 
cury’s orbit within less than a few million years 
(Figs. 2, 6). Once Mercury was removed (merged 
with Venus/Sun via inelastic collision), the sys¬ 
tem behaved very stable (no indication of chaotic 
behavior in the BS integrations). 

3.2. High-ex solutions: Symplectic vs. 

Bulirsch-Stoer results 

As mentioned above, the symplectic integra¬ 
tions of the ten high-ex solutions (smallest At = 
1/4 day) quickly failed once ex reached 0.90-0.95. 
In five out of ten runs, the deterioration of the inte¬ 


grations was immediately obvious because the rel¬ 
ative energy error, \AE/E\, rapidly grew beyond 
10“® due to the constant 1/4-day timestep, which 
was clearly too large for these runs (Fig. 6). How¬ 
ever, the symplectic \AE/E\ remained ’initially’ 
below '^10“® in the other five runs, which, if fo¬ 
cusing only on the maximum energy error, may 
not raise a red flag. Yet, large oscillations and 
jumps in the symplectic \AE/E\ occurred in these 
five runs during high-ex periods, after which the 
symplectic trajectories rapidly diverged from the 
BS trajectories (Fig. 6). (’Initially’ here refers to 
the time interval when \AE/E\ < 10“® but tra¬ 
jectories already diverge; if symplectic At is re¬ 
duced subsequently, spurious behavior may go un¬ 
noticed). Note that while \AE/E\ was still < 10“® 
after these events, \AE/E\ in the symplectic in¬ 
tegrations was typically two orders of magnitude 
larger than in the BS integrations at that point. 
In runs 0043, 0299, and 0397, ox values in the 
symplectic integrations subsequently dropped be¬ 
low those of the BS runs. In runs 1461 and 1594, 
ax shifted to larger values in the BS integrations, 
which are missed in the symplectic integrations. 
These observations reiterate the point made above 
that even though the symplectic relative energy 
error may remain below a certain threshold (say 
< 10“®), this does not guarantee accurate orbit 
integration, for instance, during close encounters 
and for highly eccentric orbits. 

Ultimately, eight of the ten symplectic high- 
ex integrations resulted in large oscillations and 
jumps in \AE/E\ to values >10“®. In addition, 
runs 0043, 0299, and 0728 resulted in rapid desta¬ 
bilization of Mercury’s orbit (ex > 0.99 and large 
shifts in ax)- Two runs (1461 and 1594) con¬ 
served \AE/E\ values to below 10“® but featured 
a rapid, suspicious decline in ex and missed the 
ax shifts seen in the BS runs (see above). Thus 
all ten symplectic high-ex integrations essentially 
failed within only ^4 Myr once ex had reached 
critical values between 0.90 and 0.95. The short 
lifetime of the highly eccentric Mercurian orbits 
in the symplectic integrations (compared to the 
much longer lifetime in the BS integrations, see 
Fig. 5) is of course mostly a result of too large a 
timestep, which was constant throughout the sym¬ 
plectic integration. However, the critical point is 
that even if the symplectic timestep was reduced 
during the integration to maintain a certain energy 
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Fig. 7.— Illustration of the stability of Earth’s orbit, (a) Earth’s orbital path in the heliocentric Cartesian system 
(shown are 15 orbits of ROOOl, separated by 100 kyr each). The area of intercept of Earth’s trajectory with the 
xz-plane (x > 0) is indicated by the purple rectangle (size not to scale), (b) Dots represent the intercepts during the 
first 1.6 Myr of R0649 plotted every 1,600 yr. (c) Intercepts during the final 2.3 Gyr of ROOOl plotted every 160 kyr. 
(d) Intercepts during a 4-Myr period of R0043 plotted every 2,000 yr. The minimum and maximum a:-values (~0.85 
and ~1.15 AU) occur during a -^200 kyr interval just before Mercury collides with Venus, (e) Poincare section of 
Earth’s trajectory in phase space: x-velocity vs. x-coordinate when Earth’s trajectory crosses the x 2 :-plane (x > 0), 
corresponding to (d). (f) Same as (e) but ^-velocity vs. ^-coordinate. Note that in all 1,600 solutions integrated here. 
Earth’s orbital path intersected the xa-plane within an area bounded by x and z-coordinates that varied at most 
by ±0.15 and ±0.05 AU from the mean. The corresponding x and z-velocities varied at most by ±2.5xl0~® and 
±9.0x10“'^ AU day“^ from the mean. 

3.3. Future evolution of Earth’s orbit 


error (At can not be changed too often though), 
symplectic integrators can easily produce spurious 
results for highly eccentric orbits and during close 
encounters, as demonstrated by the comparison 
with the BS integrations. 


Mercury’s orbital dynamics had little effect on 
Earth’s orbit as none of the 1,600 solutions showed 
a destabilization of Earth’s future orbit over the 
next 5 Gyr. Rather, Earth’s future orbit was 
highly stable. For illustration, during a typical, 
full 5-Gyr run. Earth’s orbital path typically in¬ 
tersected the xz-plane of the heliocentric Carte- 
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sian system within an area bounded by x and 
z-coordinates that varied at most by ±0.07 and 
±0.05 AU from the mean, respectively (Fig. 3.2). 
In the most extreme case found here (R0043), x 
and z varied at maximum by ±0.15 and ±0.05 AU 
(Earth’s eccentricity approaching 0.15). Impor¬ 
tantly, the maximum x, z variations in R0043 were 
actually restricted to a ~200 kyr interval just 
before Mercury collided with Venus (Fig. 3.2). 
Poincare sections of Earth’s trajectory in phase 
space showed that the corresponding x and z- 
velocities in all 1,600 solutions varied at most by 
±2.5x10“^ and ±9.0x10“"^ AU day“^ from the 
mean (Fig. 3.2). 

Also, none of the 1,600 simulations led to a 
close encounter, let alone a collision, involving the 
Earth. The minimum distance between the Earth 
and another planet (viz. Venus, dmin — 0.09 AU) 
occurred in R0728 during a 50-year period when 
Mercury collided with the Sun. This dmin does 
not qualify as a close encounter though; it is still 
^9 times larger than the Earth-Venus mutual Hill 
radius, rn — 0.01 AU (Chambers et al. 1996), or 
^1/3 of the current dmin (~30 x th)- A previous 
study suggested the possibility of a collision be¬ 
tween the Earth and Venus via transfer of angular 
momentum from the giant planets to the terres¬ 
trial planets (Laskar & Gastineau 2009). However, 
the total destabilization of the inner Solar System 
only occurred after cm Fad already crossed 0.9 
and |A£’/£’| had grown beyond 2x10“® in their 
symplectic integration (Laskar & Gastineau 2009). 
Did such disastrous trajectories for the Earth arise 
because Mercury’s periapse and close encounters 
were not adequately resolved at some point in the 
symplectic integration? 

The values of Earth’s orbital elements semi¬ 
major axis (as), eccentricity (es), and inclina¬ 
tion (is) remained within narrow bands in all 
1,600 solutions studied here (Fig. 8). For instance, 
max(a£) < 1.0005 AU, max(e£) < 0.15, and 
max(i£:) < 5.1 deg (for comparison, J2000 mean 
values are as = 1.0000, es = 0.017, is = 0.0deg). 
As the present study essentially sampled eight tril¬ 
lion annual realizations of possible Solar System 
configurations in the future (1,600 x 5 • 10®), the 
near constancy of Earth’s orbital elements a, e and 
i suggests that Earth’s orbit is dynamically highly 
stable for billions of years to come. 



Fig. 8.— Maximum values per 1-Myr bin of Earth’s 
slowly changing orbital elements. (On gyr-time scale, 
argument of periapsis, longitude of ascending node, 
and mean anomaly may be considered ’fast angles’.) 
For better visualization, only every 20th of the 1,600 
solutions are plotted plus all runs with cm > 0.8 for 
t = 2.5—5 Gyr. The 1,600 solutions used initial condi¬ 
tions that differed only by 1.75 mm in Mercury’s initial 
radial distance between every two adjacent orbits, (a) 
Maximum semi-major axis, (b) eccentricity, and (c) in¬ 
clination. All three elements remained within narrow 
bands in the 1,600 solutions: max(a£) < 1.0005 AU, 
max(e£) < 0.15, and max(i£) < 5.1 deg. The runs 
with elevated as and es are solutions with large in¬ 
creases in Mercury’s eccentricity. For example, as 
reached 1.00043 AU in R1461 during a high-ejn inter¬ 
val, which occurred ~8 Myr before Mercury collided 
with Venus. After Mercury had been merged with 
Venus, the system behaved very stable. 
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4. Conclusions 


In the present manuscript, I have reported re¬ 
sults from computationally demanding ensemble 
integrations {N = 1,600) of the Solar System’s 
full equations of motion at unprecedented accu¬ 
racy over 5 Gyr. All integrations included rela¬ 
tivistic corrections, which substantially reduce the 
probability for Mercury’s orbit to achieve large ec¬ 
centricities. The computations show that the rel¬ 
ative energy error in symplectic integrations (say 
< 10“®) is not a sufficient criterion to ensure ac¬ 
curate steps for highly eccentric orbits and during 
close encounters. The calculated odds for a large 
increase in Mercury’s eccentricity are less than 
previously estimated. Most importantly, none of 
the 1,600 solutions led to a close encounter, let 
alone a collision, involving the Earth. I conclude 
that Earth’s orbit is dynamically highly stable for 
billions of years in the future, despite the chaotic 
behavior of the Solar System. A dynamic life¬ 
time of > 5 Gyr into the future may be somewhat 
short of the Sun’s red giant phase when most inner 
planets will likely be engulfed, but clearly exceeds 
Earth’s estimated future habitability of perhaps 
another 1-3 Gyr (Schroder & Smith 2008; Rushby 
et al. 2013). 
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A. Contributions from general relativity 

Relativistic corrections (Einstein 1916) are critical as they substantially reduce the probability for Mer¬ 
cury’s orbit to achieve large eccentricities (Laskar & Gastineau 2009; Zeebe 2015). General relativity (GR) 
corrections are available in HNBody but not in mercuryS. First Post-Newtonian (IPN) corrections (Soffel 
1989; Poisson & Will 2014) were therefore implemented before using mercury6’s Bulirsch-Stoer algorithm. 

Because of the Sun’s dominant mass, GR effects are considered only between each planet and the Sun, 
not between planets. Relative to the barycenter, the IPN acceleration (denoted by ~ ) due to the Sun’s 
mass mi on planet j = 2, N may be written as (Soffel 1989; Poisson & Will 2014): 


a, = 


1 I Gmi 


J.3 

j 


9 „ 9 w N 5 , 1 , 5Gm,- 4Gmi 

-v^ - 2ui + ■ vi) + • vi) + —— + —— 

j '3 '3 






(Al) 


where a;’s and u’s are barycentric positions and velocities, the superscript h refers to heliocentric coordinates, 
and Tj = \x^\. The IPN acceleration on the Sun is: 


N 

1=2 


1 J Gmj 


Grrii 


Ui -f 2vj - 4,{vj -Vi)- -Vj) - 


2 AGrrij 5Gmi 


{4vi - 3t>j )] 


(A2) 


As mercuryS’s Bulirsch-Stoer algorithm uses heliocentric coordinates (xj = Xj — Xi), the IPN acceleration 
on planet j in heliocentric coordinates is required, given by: 


= aj — ai . 


(A3) 


The above equations were implemented in the force calculation routines of mercuryS. 

Furthermore, the energy and angular momentum correction terms in the IPN approximation are (Poisson 
& Will 2014): 


Ej = 


Lj = 




(3+ „,)(»?)>+ 

3 




i(l - 3vj){vjf + (3 -f 


X 


(A4) 

(A5) 


where rjj = mirrij/M'j and Mj = mi -|- ruj. These correction terms were added to the routine mxx_en(), 
which computes energy and angular momentum in mercuryG. Finally, the IPN Solar System barycenter in 
the present approximation is given by (Newhall et al. 1983): 


0 = m\X\ 




Gmi 
2c^ ^ 2c^ri 


(A6) 


which was used in the conversion between heliocentric and barycentric coordinates. 
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The results obtained with mercury6 and the above GR implementation may be compared to results 
obtained with HNBody (both Bulirsch-Stoer, relative accuracy 10“^®). For example, over the 21st century, 
Mercury’s average perihelion precession (only due to GR) was 0.42977” y“^ computed with HNBody and 
0.42976” y“^ computed with mercuryS and IPN corrections. In terms of Mercury’s eccentricity (ex) 
evolution, the difference in ex between HNBody and mercuryS runs (both Bulirsch-Stoer with GR correction) 
was < 10“® over 2 Myr starting at present initial conditions. 
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Table 2: Initial conditions of the eight planets and 
Pluto“ for 5-Gyr runs from DE431. Heliocentric 
positions x (AU) and velocities v (AU day“^). 



X 

y 

Z 

X 

-1.40712354144735680E-01 

Mercury 

-4.43906230277241465E-01 

-2.33474338281349329E-02 

V 

+2.11691765462179472E-02 

-7.09701275933066148E-03 

-2.52278032052283448E-03 

X 

-7.18629835259113170E-01 

Venus 

-2.25188858612526514E-02 

+4.11716131772919824E-02 

V 

+5.13955712094533914E-04 

-2.03061283748202266E-02 

-3.07198741951420558E-04 

X 

-1.68563248623229384E-01 

Earth + Moon 
+9.68761420122898564E-01 

-1.15183154209270563E-06 

V 

-1.72299715055074729E-02 

-3.01349780674632205E-03 

+2.41254068070491868E-08 

X 

+1.39036162161402177E+00 

Mars 

-2.09984400533893799E-02 

-3.46177919349353047E-02 

V 

+7.47813544105227729E-04 

+1.51863004086334515E-02 

+2.99756038504512547E-04 

X 

+4.00345668418424960E+00 

Jupiter 

+2.93535844833712467E+00 

-1.01823217020834328E-01 

V 

-4.56348056882991196E-03 

+6.44675255807273997E-03 

+7.54565159392195741E-05 

X 

+6.40855153734800886E+00 

Saturn 

+6.56804703677062207E+00 

-3.69127809402511886E-01 

V 

-4.29112154163879215E-03 

+3.89157880254167561E-03 

+1.02876894772680478E-04 

X 

+1.44305195077618524E+01 

Uranus 

-1.37356563056406209E+01 

-2.38128487167790809E-01 

V 

+2.67837949019966498E-03 

+2.67244291355153403E-03 

-2.47764637737944378E-05 

X 

+1.68107582839480649E+01 

Neptune 

-2.49926499733276124E+01 

+1.27271208982211476E-01 

V 

+2.57936917068014599E-03 

+1.77676956230748452E-03 

-9.59089132565213410E-05 

X 

-9.87686582399026491E+00 

Pluto 

-2.79580297772433077E+01 

+5.85080284687055574E+00 

V 

+3.02870206449818878E-03 

-1.53793257901232473E-03 

-7.12171623386267461E-04 


“Masses (Mercury to Pluto in solar masses): 1.66013679527193035E-07, 2.44783833966454472E-06, 3.04043264626852573E- 
06, 3.22715144505387430E-07, 9.54791938424322164E-04, 2.85885980666102893E-04, 4.36625166899970042E-05, 

5.15138902053549668E-05, 7.40740740740740710E-09. 
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